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ABSTRACT 

The standard point vortex method has recently been shown to be of high 
order of accuracy for problems on the whole plane, when using a uniform 
initial subdivision for assigning the vorticity to the points. If obstacles 
are present in the flow, this high order deteriorates to first or second- 
order. This paper introduces new vortex methods which are of arbitrary 
accuracy (under regularity assumptions) regardless of the presence of bodies 
and the uniformity of the initial subdivision. 
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1. INTRODUCTION 


There has been a growing interest recently in the theory and application 
of point vortex methods to the numerical solution of the incompressible Euler 
and Navier-Stokes equations. The impetus for the Euler case stems from the 
basic work of Dushane [6], Raid and Del Prete [7], and Hald [8], the Fourier 
analysis of Beale and Majda [1], [2], [3], and the Sobolev space approach of 
Raviart [12] and Cottet [4], A recent paper by Cottet and Gallic [5] extends 
the latter approach to linear Burger's type equations with "viscosity" 
accounted for by splitting the convection and viscous parts and using a 
Green's function for the viscous computation. A method for introducing 
viscosity into particle methods for compressible flows is given by Monaghan 
and Gingold [9]. See also [10] and [11]. Apart from the first three of these 
references, the authors all obtain high order of accuracy error estimates, 
limited mainly by the regularity of the exact solution of the continuous 
equations. Unfortunately, the possibility of obtaining this accuracy is 
dependent on the existence of expansions similar in nature to the Euler- 
MacLaurin sum formula. If, for any reason, it is not possible to assert the 
existence of such expansions, the accuracy drops to first- or second-order, 
depending on the exact details of the algorithm and which errors are being 
estimated. If general boundaries (bodies) are present in the flow field, or 
if the initial subdivision of the flow field is not uniform, the necessary 
expansions will most likely cease to exist. Then questions arise as to how 
higher-order schemes may be constructed, and more important whether it is 
worthwhile to use them in view of the extra expense which is involved. The 
purpose of the paper is to give some possible answers to these questions. 
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In Section 2, the basic equations are given, and the simplest particle 
method is defined for comparison with some higher-order schemes. These 
schemes are introduced in Section 3. There, three methods for generating 
schemes of arbitrary accuracy are provided. An appendix contains some 
technical results about solving scalar hyperbolic equations with 
distributional data. 

This paper is of an algorithmic nature and does not contain numerical 
results or precise error estimates. These will appear elsewhere. 


2. MODEL PROBLEM 

The incompressible Euler Equations in vorticity-velocity form are 

io t + (uu) x + (vw) y = 0 ^ (2.1) 

\ in IF? 

div(u,v) = 0 : curl(u,v) = w I (2.2) 

with initial condition 

m(x,y,0) = w 0 (x,y). (2.3) 

The basic ideas for constructing higher-order schemes will be shown for (2.1) 

and (2.3), with (u,v) assumed given. For these linear problems it is not 
► * 

necessary to assume that (u,v) is solenoidal. 

In this setting, we will now define the basic particle (or point vortex) 
method. Subdivide the plane into squares of side h, number the squares 
1, 2, ‘$,*. 2 . in some convenient way and define a distributional approximation 
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to u) Q (x,y) by 


2 

°W x,y ^ = l h “(vV «(*-v y-y i^ (2.4) 

i 

where (x^ , y^ ) denotes the center of the 1 th mesh square, and 6(x-x i ,y-y i > 
denotes the Dirac delta function with pole at (x^,y^). Now solve (2.1) and 
(2.2) with UQ(x,y) «• ^^^(x.y). The well known solution to the latter 
problem is the distribution 

<u h (x,y,t) = l h 2 co(x i ,y i ) S(x- X(x i> y j .;t), y - Y(x i ,y ; .;t)) (2.5) 

i 

where X(x^,y^,t) denotes the solution of the characteristic equation 

dX/dt = u(X,Y,t) X(0) = x ± 

and correspondingly for Y. 

No use is made of the uniformity of the mesh in deriving (2.5). For a 
2 

nonuniform mesh, h in (2.5) is the area of the appropriate mesh square. In 
the error formulas below, h denotes the largest mesh length. 

It is immediately clear from this definition that the particle 
approximation is non-dissipative , in the sense that no artificial viscosity is 
introduced because after the discretization of the initial condition is made 
(2.1) is solved exactly. In practice some ODE solver must be used to compute 
the trajectories, but in theory its error can be made arbitrarily small. This 
principle, of solving the exact equation with approximte data, seems to be 
common to particle methods generally and distinguishes them from finite 
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difference and finite element methods. The latter, at least, solves an 
approximate equation with exact data. 

A rigorous error analysis of the method just defined can be found in 
[12]. This analysis is too complicated to reproduce here. Nevertheless, we‘ 
need some simple guide to compare the accuracy of various schemes. It seems 
reasonable to look at the difference Wq - against a test function as a 

measure of "truncation error" since it is the only error made. Thus we 
define, for a given method of approximation and a given function with 

compact support ft (where area (ft) = 1 say) 

T h (<|>) = // (“o " “0h^ dxdy * (2.6) 

Here, the integration is performed over IT. The restriction that has 

compact support is a matter of convenience rather than necessity and could be 
replaced by sufficiently rapid decay at large distances from the origin. 

As an example, consider (2.4), Then we find 

T h (<i>) = JJ “o ^ dxdy “ I h 2 (2.7) 

i 

This shows that a midpoint rule numerical integration is being used to 
approximate the integral, and under smoothness conditions it follows that as 
h + 0 

x h U) = 0(h 2 ). 

Clearly, higher-order integration formulas can be compared with each other on 

this basis. For a 2 x 2 product Gauss rule in each element, for example, we 
A 

have = 0(h ) . 
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Next, recall the important fact that in the nonlinear case it is 

necessary to compute the velocity field at each timestep by solving (2.2). 

Assume that this is to be done using the Green's function. Let W denote the 

number of arithmetical operations required to compute the velocity field at 

2 

each particle position. If there are N particles, then W s CN /2, for 
some constant C. Below, we will use W as a standard unit of work to 
compare various new algorithms. For the Gauss case therefore we have a work 
count of 16W. From this we see that use of a higher-order rule does not 
necessarily assure a greater computational efficiency for typical values of 
h. In the next section, methods for obtaining high-order accuracy without 
such a large increase in the cost of the computation are defined. 


3. HIGHER ORDER METHODS 

The preceding remarks suggest that increasing the order of accuracy by 
adding more integration nodes may not be a good idea. It is natural to try to 
do the same thing by increasing the amount of information associated with each 
node. Specifically, in this section we shall associate with (x^,y^), m 1 "' 1 
order distributions of the form 

M i (x,y) = l w ia D° 6(x-x i , y- y ± ). (3.1) 
| a j <m 

In (3.1), which generalizes the simple 6 functions in (2.4), a denotes a 
multi-index, and (x^,y^) e l£. Choice of the weights w^ a and the nodes 
(xi,yi) can be made in many ways. We shall give three methods in this 


section 
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Method 1 (Direct Integration) : 

In this method, (x^,y^) are the corners of the elements, each of which 
has associated with it an expansion of the form (3.1). The weights in the 
expansion are chosen so that when is substituted into (2.6), the second 

term gives a rule for integration of the function (wq <f>), involving its 
values along with those of its derivatives through order m at the nodes. We 
shall consider the cases m = 0 and m = 1 in more detail. 

Let m = 0. A rule for a square of side h with corners at P, Q, R, S 
which is exact for bilinear functions is 

ff f dxdy = (h 2 /4) (f(P) + f(Q) + f (R) + f(S)). (3.2) 

2 

Using this as a composite rule implies the choice w ^qq = ^ s0 

that we define 

2 

M ± (x,y) = h m(x i ,y i ) 6(x-x ± , y-y ± ). (3.3) 

Since this gives a rule which is locally exact for linear functions but not 
for all quadratics its accuracy is O(h^) in the sense of (2.6) while the 
work is 1W. This is essentially no different from the mid-point i^ile. In 
fact this rule is clearly analogous to the trapezoidal rule. 

For a quadrilateral mesh, a bilinear mapping can be used to map the 
quadrilaterals onto a standard square in which (3.1) can be used. In some 
circumstances it may be desirable to use a triangular mesh instead of the 
quadrilateral one. An O(h^) rule for triangles analogous to (3.1) can then 
be used, avoiding the need to map the domains. 
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Now let m = 1. Analogous to (3.2) we have the formula 


// f dxdy » A( f (P) + f (Q) + f(R) + f ( S) ) 

+ B(-f x (P) + f x (Q) + f x (R) - f x (S)) (3.4) 

+ C(-f y (P) - f y (Q) + f y (R) + f y ( S) ) 

where A = h^/4, B = C = h^/24, and P, Q, R, S denote the corners of the 
square -h/2 x, y < h/2 labelled counterclockwise starting from the top 
right. Analogous to (3.3) there is the expression 


M i (x,y) = l w ia D a 6(x-x 1 , y-y ± ). 

I«I<1 


(3.5) 


In (3.5), the coefficients are computed from the composite rule based on 
(3.4). For the uniform square mesh we are using for illustration, the weights 
are 


w 


iOO 


= A'« 0 (x 1 ,y i ) + B'm 0x (x i ,y.) + C' U (x.,y.) 


W ilO = 


w ioi = -c'^o^i.yi). 


(3.4) is exact for cubic polynomials. It follows that this method is accurate 
in the sense of (2.6) to O(h^). To compute work units for this scheme, we 
observe that although there are only = N particles there is some extra work 
associated with computation of derivatives of the velocity kernel. It turns 
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out that for this scheme the work units are < 2 W, a satisfactory figure. 
There is also some additional work required for computing the coefficients of 
the derivatives in (3.1). This amounts to having to integrate two more 
systems each of two odes, in addition to the characteristic odes (see 
appendix) . 

As in the previous case, rather than use a quadrilateral mesh it might 
sometimes be better to use a triangular one. 

For a square mesh, the m = 1 scheme just discussed has an interesting 
property in the uniform case. This is the following: due to cancellations, 
the composite rule has weights of zero attached to the derivative unknowns at 
interior vertices. Hence the higher accuracy is achieved by corrections at 
the boundary. But this implies the use of a Euler-Maclaurin type expansion. 
Thus, if u(f> has s continuous derivatives in I? and compact support, by 
using nodal derivatives up to this order we can get accuracy 0(h ) merely 
by using the m = 0 scheme, since this is what the composite scheme reduces 
to on a uniform mesh in that case. This is another way to look at the results 
of [1] - [3]. 


Method 2 (Finite Element Approach) : 

The approach here uses a nodal finite element basis in the following 
way: let {ij; } ! a l £. m > i = be the standard nodal basis functions 
associated with the i th node (x^,y^) of a triangulation of the plane with 
maximum edge length h. These functions satisfy conditions of the form 


tJ». (x. 
r ia v 2 


•V ' V 



-9- 


ctB 

where is a Kronecker delta. Then we define w. as 

ij ia 

w ±a = (-1)^// ^ ia (x,y) w Q (x,y)dxdy (3.6) 

where the integration is over the whole plane. We now have 

If w ot/ x,y ^ P( x »y> dxd y = // I l w ia D a 5(x-x i} y-y i ) 

* |a|<m 

x $(x,y)dxdy, V ({> € C™ (®?) 

(3.7) 

■ I I (-i)^ w ia o a ♦<* 1 iy 1 > 

1 |a|<m 

= fl to 0 (x,y) <|> h (x,y)dxdy 

where is the finite element interpolant of on the given 

triangulation. Equation (2.6) then becomes 

T h (<f>) = // “o^ ~ 4> h ) dxd y* (3.8) 

Since the error |4> - <f> h | is formally 0(h r+ ^) where r is the degree of 
the highest order full polynomial space used, we can say here that is of 

this order. 

This type of scheme differs from direct integration schemes in that no 
approximation of Wq is made. The test function only (often a convolution 
kernel in practice) is approximated and the result is integrated exactly. 
Because of this property, the rigorous error estimates for these methods 
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require minimal regularity on u)q unlike the direct integration methods where 
to achieve high accuracy requires to have several smooth derivatives 

throughout I?. The 0(h r+ *) estimate is in fact valid even if we know 

only i 0 q£L*(!£). If has extra regularity it can be exploited to get 

higher accuracy by going to negative norm estimates of the finite element 
error. Smoothness of <j> , however, is certainly required. 

Two examples analogous to those considered above are the case of 

continuous linear elements on triangles, for which we can expect O(h^) 
accuracy with 1W work units, and full cubics - defined in terms of 

derivative unknowns at vertices, and function values at vertices and centroid 
for which the work will be somewhat larger than the values used so far (about 
10 j W units) . 

In general, the full range of finite element spaces is available for use. 

Method 3 (Taylor/Moment Expansions) : 

Here we begin by subdividing the plane into arbitrary elements with mid- 
side nodes and arbitrary element shapes allowed in principle. Next, we define 

II a i a o 

a l ! a 2 ‘ w ict = a //(x-x^,) (y~ y i ) a) Q (x,y)dxdy (3.9) 

in which (x^,y^) is an arbitrary point within the i 1 "^ 1 element, and the 
integration is over the ith element. The w^ are proportional to the 
moments of u)q restricted to the i 1 "' 1 element, about (x^,y^). It follows as 
above, that 

// oi 0h (x,y) <j>(x,y)dxdy = // u) Q (x,y) <j> ^ (x,y)dxdy 


(3.10) 
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where <jJ m ^(x,y) is the piecewise polynomial function, in general 
discontinuous, equal in the i^ element to the Taylor expansion of 4>(x,y) 
through m cn order terms, about the point (x^,y^). In this sense the local 
moment expansion defined by (3.1) and (3.9) "dualizes" into the local Taylpr 
expansion about (x^,y^). 

To get the accuracy of this scheme, we substitute into (2.6) to find that 

T h u) = II _ < t ,tm])dxd y 

so that denoting by h the largest linear dimension of the elements, we 

obtain accuracy 0(h m+ ^). 

The moments method also needs only minimal regularity on idq for full 
accuracy to be obtained. In practice, if m = 1 the point (x^.y^) should 
be chosen to be the center mass of u)q because then = 0 for |a| = 1, 

so we get second-order accuracy for the same work as with the lowest-order 
scheme. Using quadrilaterals for elements, with N vertices there are 

approximately N elements and so N particles. For 0(h) accuracy the 
interaction work count is 5W, and for 0(h <!f ) is 8W. 


4. FURTHER REMARKS 

There should be no difficulty in extending the ideas of Section 3 to 
three-dimensional particle methods of the kind suggested in [1] - [3] and 
[ 12 ]. 

Rigorous analysis using the Sobolev space setting has been carried out 
for both the finite element and moment expansion methods. 
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So far an insufficient amount of computation has been done to verify the 
error estimates and decide about the efficiency of the various methods. 
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APPENDIX 


A framework for finding distributional solutions of (2.1) with initial 
condition = D a 6(x-Xq, y-yp) |a| _< m can be obtained starting from the 

following considerations. Let X(xg,yg;t) and Y(xg,yQ;t) denote the 
characteristic curves of the equation (2.1); here, t parameterizes the curve 
and the generic point (xg,yg) denotes its origin at time t = 0. X and Y 
are computed by solving the ordinary differential equations 

dX/dt = u(X,Y,t) dY/dt = v(X,Y,t) 

X(0) - x Q Y(0) = y Q . 

At time t, let J(xq, yg; t) denote the Jacobian of the flow map 
$ : (xg,yg) + (X,Y). The (nonlinear) case of most interest from the fluids 
viewpoint has u x + Vy = 0, in which case J(xg,yg;t) = 1. We can obtain a 
formal analytical solution to (2.1) and (2.3) by writing the equation in terms 
of the material derivative as du/dt = 0, integrating this equation over an 
arbitrary domain moving with the velocity field (u,v), say fi(t), and then 
using the transport theorem to write 

d/dt // w(X,Y,t)dXdY = 0, 
n(t) 

from which it follows immediately that 

// u)(X,Y,t)dXdY = // a) n (x,y)dxdy. 

n(t) D(0) 
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Changing the variables on the right-hand side to X and Y respectively and 
recalling the arbitrariness of ft(t) now gives 

u(X,Y ; t) = u 0 (x(X,Y,t), y(X,Y,t))j _1 (X,Y;t) (A.l) 

where (x(X,Y,t), y(X,Y,t)) is by inverting the equations X = X(x,y:t), Y = 
Y(x,y;t). The existence of a unique solution to these equations follows from 
ode theory provided u and v are smooth. Reversing the steps, it follows 
that (A.l) satisfies (2.1) given the required regularity of u, v, and 

Let <j i £ C m (]F?); multiplying (A.l) by <$>, integrating and changing the 
variables on the right to x and y we have 

// w(X,Y,t) 4> ( X , Y ) dXdY = // w^x.y) 4> (X(x, y ; t ) , Y(x,y,t))dxdy, (A. 2) 

or alternatively 

<u>, <p = <m 0 , <j> o(X,Y)> (A. 3) 

where o denotes composition. If X(*,*, t) and Y(*,*, t), 

Y( • , • , t) £W m+1, ”(l^) (or £ C (m) (3B?)), V 0 < t < T, then the right-side of 
(A. 3) makes sense even if Wq + = D a <5(x-Xq, y-yQ)|a| < m. Thus a 

distribution a) is defined on (f^ m ^(I^) by (A. 3). Therefore, we can pose 
the problem of finding satisfying 


<10^, *}>o(X,Y)> 


¥ <j> £ (l£). 


(A.4) 
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A solution o)^ to (A. 4) is given by 

u. (X,Y) = D a 6(X - X(x,y;t), Y - Y(x,y;t)) , (A.5) 

n x-x Q , y-y Q 

the purely formal differentiations being performed w.r.t. x and y. Proof 
that (A.5) satisfies (A. 4) is by direct computation. 

If | a | = 0 we recover the solution given in Section 2. Consider the 
case with | ex | = 1. Equation (A.5) gives 

“hio ■ M x - x 0' Y - Y o> MW 1 ) + M x -V Y - Y o) Y x ( X 0' Y 0- £ ) 

(A. 6) 

“hOl ■ S x( X -V Y ‘ Y o) X y^ x 0’ Y 0’ + M X - X 0' Y - Y o) V W'l 

using the abbreviation Xq for X(xQ,yQ;t) and similarly Yq. If the 
initial condition is • 

“ho = a io 6 x( x " x o> y_y o) + a oi M x - x o> y ' y o)’ 

then the solution to (A. 4) of the required form as given by (A. 6) is 

u h = a io (t) 6 x^ x “ x o> Y " Y o) + a oi (t) y_ V 

where 

a l0 (t) = a 10 MW') + a 01 V 

(A._7) 

a 0i(t) = a 1Q Y x (x 0 ,y 0 ,t) + a Q1 Y y (x 0 ,y Q ,t). 
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Letting M denote the matrix 



differentiation of the characteristic equations shows that 


dM/dt = V(u,v)M 


and the initial condition for this system is M(0) = 1, the identity matrix. 
It will be necessary to solve this and analogous systems for the higher-order 
cases in order to compute the numerical approximations. Having solved it, 
a^Q(t) and agj(t) are given by (A. 7). 
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